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1 - INTRODU CTION 

In recent years a multitude of techniques has been developed for 
generating computational grids required in the finite difference or 
finite element solutions of partial differential equations on arbitrary 
regions. The importance of the choice of the grid is well known. A 
poorly chosen grid may cause results to be erroneous or may fail to 
reveal critical aspects of the true solution. Some considerations 
that are involved in grid selection can be noted from the papers of 
Blottner and Roache [1], Crowder and Dalton [2], and Kalnay de Rivas 
[3]. While these papers discuss error for one-dimensional problems, 
few results exist for higher dimensions. This report will examine the 
errors in approximating the derivatives of a function by traditional 
central differences at grid points of a curvilinear coordinate system. 

The implications concerning the accuracy of the numerical solution of 
a partial differential equation will be explained by considering several 
numerical examples. Although this study only considers the two-dimensional 
case, the techniques and implications are equally valid for three- 
dimensional grids. 

An interesting feature of the error analysis in this report is its 
simplicity. Most of the results follow by merely working with the trunca- 
tion terms of some power series expansion. It is noted that these series 
expansions also give rise to higher-order difference approximations 
which can significantly reduce error when the grid spacing changes rapidly, 
as might be the case in problems with shock waves or thin boundary layers. 
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When transforming a partial differential equation from rectangular 
to curvilinear coordinates, the derivatives of the functions defining 
the transformation must be evaluated. If the relation between rectangular 
and curvilinear variables is given by a simple analytic expression, the 
transformation derivatives may be computed either analytically or 
numerically. Truncation errors in both cases are considered for 
comparison . 

One objective of this work is to provide tools to examine a grid, 
together with a computed solution, and predict possible inaccuracies 
due to the grid. The grid may thus be redefined to give a better 
solution. Directions for future work could be an extension to higher 
dimensions of the one-dimensional grid optimization technique of Pierson 
and Kutler [4], 

This report also discusses the control of coordinate line spacing 
through functions incorporated in the elliptic generating system for 
the curvilinear coordinates. Attraction of coordinate lines to other 
coordinate lines and also attraction to fixed lines in physical space 
are covered. Appropriate forms of the control functions required to 
produce desired spacing distributions are derived. Finally a procedure 
for distribution of points around a boundary curve according to local 
boundary curvature is given. In addition a few examples of recent 
generation of coordinate systems are given. 


II . TRUNCATI ON ERROR AN ALYSIS 


Suppose a curvilinear coordinate system is generated by trans- 
forming an arbitrary physical region of the xy-plane onto a rectan- 
gular computational region of the £n-plane. The relationship between 
partial derivatives of a function f with respect to physical and 
computational variables is well-known. It will be included here for 
later comparison with approximations derived from series expansions. 
Only first and second order derivatives will be considered: 
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The derivatives with respect to p can be obtained by replacing £ with 
p in the first two equations of (1). 

Although this change of variables formulation can be easily used 
in deriving difference approximations for derivatives with respect to x and 
nothing can be said about truncation error. An error analysis can, however, 
be based on Taylor series expansion of function values at neighboring points 
about a single point in the physical region. In order to distinguish 
between derivatives and differences in the following, the dif f erent ial 
notation is used for derivatives while subscripts denote the usual 
second order central difference expressions . The following approximations 
for the central differences are valid when all series are truncated after 
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second derivative terms. A unit mesh width in the £n~plane is assumed 
without loss of generality. 
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Together with the corresponding two equations for f n and fqq, this 

constitutes a system of five simultaneous equations which can be solved 

to produce difference expressions of two first and three second derivatives 

of f with respect to x and y. Assuming the third order derivatives of 

3 

f are bounded, the truncation error in the above expressions is 0(h ), 

where h is some measure of the local mesh spacing. Consequently, when 

(2) is solved for the difference approximations of the physical derivatives 

2 

of f, the truncation error is 0(h ) for first derivatives and 0(h) for 
second order derivatives. In contrast, solving the system (1) with the 
f, x, and y derivatives replaced by differences (and including the 
corresponding equations for f n and f rjr>) simultaneously to produce expres- 
sions for the five physical derivatives of f gives rise to 0(h) and 0(1) 
truncation errors for the first and second order derivatives. 

In both cases it has been assumed that the coefficient matrices on 
the right hand sides of (1) and (2), i.e., the coordinate derivatives, 
including the omitted n differences, are well conditioned. Ill conditioned 
matrices which may result from extremely skewed coordinate lines could 
cause further deterioration in accuracy. Higher order accuracy can be 


obtained using (1) if second order coordinate differences are assumed 
2 

to be 0(h ). This effectively limits the rate of change in coordinate 
line spacing and the curvature of coordinate lines, however. No 
simple relation between the coefficients of the second order derivatives 
in the last equation of (1) and (2) was found except for the fact that 
they would be equal if the differences in (2) were replaced by deriva- 
tives. 

The variation in numerical solutions using (1) and (2) is illustrated 
in the solution of Laplace* s equation. The function 

u(x, y) = x(l + l/(x ? + y 2 )) 

2 2 

satisfies Laplace* s equation for x + y >1 and has a vanishing normal 

derivative on the boundary. This boundary value problem was solved 

2 2 

numerically on 1 < x + y < 100. A grid was selected with 39 radial 
coordinate lines and 49 circular coordinate lines. The first 23 
circular coordinate lines were uniformly spaced after which the spacing 
was increased by a factor of 5. The difference between the exact and 
numerical solution is indicated in Figure 1 for difference equations 
derived from (1) and (2). The effect of the sudden change in coordinate 
line spacing was clearly less severe when using difference expressions 
from the higher order series expansion. 

A similar error analysis can be carried out where the derivatives 
of x and y with respect to £ and n are computed analytically rather than 
approximated by differences. In this case a series expansion in the 
£n-plane is required, followed by substitution of expressions for the 
higher-order £ and n derivatives in terms of the x and y derivatives 
(see Ref. 5 for complete detail). Retaining physical derivatives of f 
through second order, as in (2), the following approximations are 
generated. The second derivative approximations, f^ and f^rp are ver Y 
lengthy and only the first and second order derivatives of x and y are 
included here, the complete expressions being given in Ref. 5. The 
first derivative approximation includes third order derivatives: 
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Considerable similarity exists between the approximations in (2) 

and (3) and corresponding statements can be made about the effects of 

the coordinate system on truncation error. For example, it can be noted 

that for the first derivative approximations to be second order accurate, 

2 

the second and third order derivatives of x and y must be 0(h ) and 

3 

0(h ), respectively. Due to the additional restriction on the third 
order derivatives, it is not difficult to find examples where solutions 
of (1) with numerically computed derivatives of x and y are much more 
accurate than solutions using the analytical expressions for these 
derivatives. 

With reasonable care in the selection of the grid any of the above 
difference formulations will give equally good results. For example, 
consider the grid for the region about a Joukowski airfoil depicted in 
Figure 2. This grid was constructed by the conformal mapping of an 
annular region with uniformly spaced circular coordinates. As in the 
above example, Laplace’s equation is solved with vanishing normal 
derivative imposed on the airfoil. The solution is the velocity potential 
for flow about the airfoil at zero angle of attack. Table 1 indicates the 

v. 
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the difference between the computed solution and the exact solution 
on the surface of the airfoil where the error was greatest. 


Table 1. Comparison of Difference Formulations 


Differencing Method 

Max Error 

RMS Error 

Taylor Series (2) 

.03123 

.00864 

Analytic (1) 

.02216 

.01256 

Numerical (1) 

.02411 

.00795 


For this example there is clearly no advantage in using the difference 
expression from the series expansion in (2) over using (1) with the 
derivatives of x and y computed either analytically or numerically. 

There is another aspect to the question of the use of analytically 
calculated coordinate derivatives, as opposed to numerical difference 
representatives, when fully conservative difference formulations are 
used. In that case the formulation will not be fully conservative with 
the analytical expression in the sense that a uniform solution on the 
field will not be strictly preserved. This can lead to instability if 
the differences of the coordinate derivatives are large. 

Thus far only problems of error which deal directly with the coor- 
dinate system have been considered. This source of error can be controlled 
by limiting the higher order differences of derivatives of x and y. A 
more serious problem in numerical computations is the error in the approxi- 
mate solution which results from large higher order derivatives of f. 

In transforming from physical to computational variables, the derivatives 
of f with respect to £ and r| are replaced by differences regardless of 
whether derivatives or differences are used for x and y. The truncation 
error in approximating the computational derivatives of f can be minimized 
to some degree by a properly chosen grid. However, there are limitations 
in the grid choice since, as we have previously observed, a highly dis- 
torted grid also contributes to large truncation errors in the approxima- 
tion of the physical derivatives of f. To analyze the total truncation 
error due to solution and grid, it is convenient to introduce matrix 
notation. 

Suppose the derivatives in the physical and computational planes 
are related by (1) . This relation can be written 


(M 


where 


A * AD 







Difference expressions for D are generated by replacing the elements 
of A (and possibly A) by the appropriate difference approximations. If 
the truncation term is retained, the equation (4) becomes 

6 + e = AD (5) 

where 6 is the vector containing the difference approximations and 
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Solving for D in (5) we have 

D = A -1 <5 + A _1 e 

Now e is unknown but can be estimated using differences of f. 
Although such numerical differentiation does not tend to be very 
accurate when applied to an approximate solution of a partial differ- 
ential equation, the value of A has been used successfully to 
distinguish regions of high error from regions of low error. This can 
be illustrated by returning to the numerical solution of potential flow 
about the Joukowski airfoil in Figure 2. The comparison of truncation 
error with error in the solution is indicated in Figure 3 for grid 
points beginning near the trailing edge and ending near the leading 
edge of the airfoil. The grid points were chosen to lie on the second 
coordinate line from the airfoil surface so that no extrapolation was 
needed to estimate the elements of e. 

Each factor in the truncation error estimate can be analyzed 
independently. The factor A 1 deals only with the grid coordinates, 
while c involves only the solution of the partial differential equation 
In the above example consideration of e alone would seriously under- 
estimate the order of accuracy near the leading and trailing edge since 
the distortion in the coordinate system would not be taken into account 
The influence of the factor A ^ can be analyzed by examining the condi- 
tion of the matrix A. An ill-conditioned matrix not only magnifies the 
effect of the truncation terms In e but also the effect of deleting the 
additional terms which appeared in the series expansions (2) and (3). 

We will now consider a case where an extremely ill-conditioned 
matrix is encountered. The Navier-Stokes equations in stream function- 
vorticity formula t ion were solved numerically for viscous flow about a 
circular cylinder. The data in Table 2 illustrates the growth in the 
condition number of A as the circular coordinate lines are concentrated 
near the cylinder to resolve the boundary layer. Only the Laplacian 
of vorticity was included in the truncation error computation since 
this truncation term clearly dominated the remaining truncation terms 
in the equations. As n increases, the dominating factors in the trunca- 
tion term shifts from the elements of e to the elements of A \ An 
examination of vorticity values revealed a clear deterioration of the 
numerical solution for n * 4. 


Table 2. Maximum truncation error for V 2 w and 
condition number of A. Circular coordinate lines 
r - 1 -f 9(1 - exp(nn/48)) / (1 - exp(n)), 
r) = 0, 1, , 48. Reynolds no. = 5. 


n 

Max Truncation 

Max Cond^ (A) 

1 

.1079 

45 

2 

.0482 

78 

3 

.0891 

259 

4 

3.0918 

1017 


For later reference, we have from (2) for the one-dimensional case 
that the simple two-point central difference expression for the first 
derivative, f^/x^, has a truncation error term given by 






which acts as a numerical diffusion. This effect was pointed out earlier 
in Ref. 1. 



III. COORDINATE SYSTEM CONTROL 


A. O riginal Generating System 

In the formulation of boundary-fitted coordinate systems generated 
from elliptic systems as given in Ref. 6 the curvilinear coordinates (£, ri) 
were determined as the solution of the system 


V 2 ^ = P(£,n) (7a) 

v 2 n = QU,n) (7b) 


which in the transformed plane becomes (from here on, s ubscripts indicate 
derivatives ) 


ax 




2$X__ + yx 


nn 


= -J (Px + 




(8a) 


- 2 % n 


+ yy 


nn 


= -J (Py^ + Qy n ) 


(8b) 


with 


a 


2 

x 

n 



(9a) 


e i 


x^x 

4 n 


Vn 


(9b) 


47 


(9c) 
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B . Attraction to Coordinate Lines 

Here the functions P and Q are to be chosen to control the coordinate 
line spacing. In Ref. 6 those control functions were taken as sums of 
decaying exponentials of the form 


P = -z a.sgn (£ - £ )exp(-c.|s “ £,|) 
i=l 1 1 11 
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Here the a., b., c., and d, of the Q functions are not necessarily the 
i i i i y 

same as those in the P function. 

In the P function the effect of the amplitude is to attract £- 
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coordinate lines toward the £^-line, while the effect of the amplitude 
tK is to attract £-lines toward the single point n^) ■ Note that 

this attraction to a point is actually attraction of £-lines to a point 
on another £-line, and as such acts normal to the £-line through the 
point. There is no attraction of p-lines to this point via the P 
function. In each case the range of the attraction effect is determined 
by the decay factors, c^ and d^ . With the inclusion of the sign changing 
function, the attraction occurs on both sides of the £-line, or the 
(£f, r\ ) point, as the case may be. Without this function, attraction 
occurs only on the side toward increasing £, with repulsion occur ing 
on the other side. 

A negative amplitude simply reverses all of the above-described 
effects, i.e., attraction becomes repulsion and vice versa. The effect 
of the Q function on q-lines follows analagously. A number of examples 
of this type of coordinate line control have been given in Ref. 6. 

In the case of a boundary that is an n-line, positive amplitudes 
in the Q function will cause n-lines off the boundary to move closer to 
the boundary, assuming that n increases off the boundary. The effect 
of the P function will be to alter the angle at which the £-lines inter- 
sect the boundary, since the points on the boundary are fixed, with 
the £-lines tending to lean in the direction of decreasing If the 
boundary is such that n decreases off the boundary then the amplitudes 
in the Q function must be negative to achieve attraction to the 
boundary. In any case, the amplitudes a ^ cause the effects to occur 
all along the boundary, while the effects of the amplitudes b_^ occur 
only near selected points on the boundary. 

If the attraction line and/or the attraction points are In the 
field, rather than on a boundary, then the attraction is not to a fixed 
line or point in space, since the attraction line or points are them- 
selves solutions of the system of equations, the functions P and Q 
being functions of the variables £ and p. It is, of course, also 
possible to take these control functions as functions of x and y, instead 
of £ and p, and achieve attraction to fixed lines and/or points in the 
physical field. This case becomes somewhat more complicated, since it 
must be ensured that coordinate lines are not attracted parallel to 
themselves, and its discussion follows in a later section. 
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C . Control Functions for Certain Spacings 

For certain simple geometries it is possible to integrate (8) 
analytically for appropriately selected forms of the control functions, 
and thus to determine the control functions required to produce a certain 
line spacing. In this regard consider the case of two concentric circu- 
lar boundaries of radii r^ and r 2 > with r 2 > rj. 

With n = 1 on the inner boundary, n * J on the outer boundary, and 
£ varying monotonically from 1 to I around these boundaries, a solution 
of (8) can be given in the form 


x = r(n) cos [2 tt 


<&» 


(11a) 


y * r(n) sin [2tt (j^) 3 


(lib) 


Substitution of these expressions into the equations of (8) with P(C,n) = 0 
yields 
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ii 


r' 



r ' 2 Q = 0 


( 12 ) 


This can be made a perfect differential by taking the control function 
Q to be of the form (following the direction of Ref. 7) 
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f"(n) 1 

f'(n) r , 2 


where the minus sign has been introduced merely for convenience. Since 
1 Y 

— =■ is equal to -hr for the solution given by (11), this form of Q suggests 
r' 1 J 

taking Q to be of the form 
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Substitution of (13) into (12) yields 


f" 

f - 0 


( 14 ) 


which can be integrated twice to yield 

c x f (n) 

r(n) = C 2 e 


The constants of integration may be evaluated from the boundary condi- 
tions, r(l) = r^, r(J) = ^ » so that 

f f(n) - fOl l 

r(n) = r 1 {(^/r^ f(J) " f(1) } (15) 


This equation may then be solved for f(n) to yield 


f(n) - f(l) 
f(J) - f(l) 
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If the distance from the body to the Nth inline is specified to be 
the following equation must be satisfied: 


f(N) - f(l) 
f(J) - f(l) 


ln[- 


r 2 

In [ — — ] 
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(17) 


It should be noted that the form of f(n) is still arbitrary, subject to 
(17). 



Alternatively, to set r* at the inner boundary, n = 1 , we have, 
upon differentiation of (15) with respect to n and subsequent evaluation 
at n - 1, that f(n) must satisfy 


iim. 

fa) _ r i 
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(18) 


The two derivatives appearing in the truncation error of first 
derivatives, as given in last equation in section II, are, from repeated 
differentiation of (15), 
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Thus if a function f(n) with a free parameter is selected, (17) 
may be used to determine the parameter in the function such that the 
Nth n-line lies at a specified distance, r^ - r^, from the inner 
boundary. Alternatively, the free parameter may be determined by (18) 
such that the spacing at the boundary is set by specification of r T 
there. The derivatives in the truncation error terms may then be 
calculated from (19). With the function f(n) determined, the control 
function Q is then given by (13). 

For example, with the function (Ref. 7) 

f(n) = nK n-1 


where K is a free parameter, we have, by (17), that K must be the 
solution of the nonlinear equation 
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to set the Nth r\-line at . Alternatively, the value of K required 
to set r 1 to a specified value at the inner boundary is determined by 
(18) as the solution of the nonlinear equation 


1 + InK 
JK J_1 - 1 


r ’ (1) 
r. 


r 2 

ln(— — ) 
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( 21 ) 


For this function, the derivatives appearing in the truncation error 
term, (19) , are given by 


f ’ (n) = (1 + nlnK)K n_1 


(22a) 


f"(n) = (2 + nlnK)(lnK)K n 1 
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The control function Q is given by (13) as 


Q - 
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j2 V l + nlnK 
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It can be shown by consideration of the ratios of successive 
derivatives that the higher derivatives of this function are progres- 
sively decreasing if K is in the range 


0 < InK < - 1) 


(24) 


Since the left side of (21) is a decreasing function of K for positive 
K, the smallest value of the spacing at the boundary, r 1 (1), that can 
be achieved while maintaining progressively decreasing higher derivatives 
of f (n) occurs with K at the upper limit of the inequality (24), viz 


r ? y(l + /5) 
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It is not reasonable to use smaller values of r r (l) since the progressively 
increasing higher derivatives of f(n) will result in significant trunca- 
tion error introduced by the coordinate system. 

Another choice of f(n) might be 


f(n) = sinh [K(n - 1)] 


(26) 


for which, from (17), the Nth line occurs at r^ 
solution of the equation 


sinh [K(N - 1)| 
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for K given by the 


or the spacing at the inner boundary is r ! (l) for K given by (18): 
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The first two derivatives and the control function are given by 


f 1 (n) - Kcosh[K(n - 1)] (28a) 
f"(n) = K 2 sinh [ K(n - D] (28b) 
Q = - ^2 Ktanh (K(n - 1)] (29) 


In this case progressively decreasing higher derivatives occur for K in 
the range 0 < K < 1, so that the smallest practical spacing at the 
inner boundary is 


r 1 (1) 


r 2 
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The control function for the spacing distribution of Roberts, 
Ref. 8, can be determined in the same manner as follows. With the 
notation of Ref. 8 adjusted so that the boundaries occur as used 
above, we have 


r(n) = r 2 
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c(n) - 1 
G(n) + 1 



(30) 


wi th 


b + r (j-~) 
G(n) = (r --) J 

b - r ~ 


with b a free parameter. 

Although the form of f(n) could be extracted by substitution of (30) 
into (16), it is simpler to determine the parameter b from either r(N) = r^, 
or for a specified value of r ! (l). The derivatives are 
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The control function Q is then given by (13) and (14) as 


(31a) 


(31b) 


Q = 



(32) 
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with r, r ? , and r M to be substituted from (30) and (31). 

Finally, another type of function is a patched function using 
different functions near and away from the inner boundary to achieve a 
group of closely spaced lines near the inner boundary with fairly 
rapid expansion outside this inner group. This is done as follows: 

Let the spacing of the inner group be such that the same change in 
velocity would occur between each two lines for a velocity distribution 
given by u(r). To do this, invert the velocity function such that 
r = r(u), and then take 

u(n) = n~^T u n i < n ^ n 


when u^ is the velocity at the edge of the inner group of lines. Then 


r(n) = r (u 


n - l 
N - 1 *V 


l i n < N 


From this function all the derivatives and the control function may be 
calculated, the latter being determined by (13) and (14). 

Now outside the inner group of lines, i.e., for N < n < J, let 
r(n) be a quartic polynomial: 

r(n) = r' N (n - N) + yr" N (n - N) 2 + |r , " N (n - N) 3 
+ a (n - N) 4 + r N N < n £ J 

where the three derivatives at n “ N are determined from the derivatives 
of the inner function, all being evaluated at n = N. The final parameter, 
a, is determined such that r(J) = r^- Thus 

(r j " V “ r V J ~ N) ■ i r Y J ■ N)2 _ l r ” ,(J ■ N)3 
a , — — 

(J - N) 4 
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The outer control function is then determined from (13) and (14). This 
composite control function has only one continuous derivative, and thus 
could possibly lead to truncation error introduced by the coordinate system. 



It is, of course, also possible to integrate the coordinate 
equation (8) for the one-dimensional case. In that case the control 
function Q is given by 


Q 


Y r M Y_ f "(n) 
j 2 r' - ' f’(n) 


( 33 ) 


and 


r(n) 


r 2 + (r 2 


" r l>f 


f(n) - f(l) 1 
f(J) - f(l) J 


(34) 


All the other steps follow in analogy with the two-dimensional case. 

Now, although the two-dimensional case given above applies only to 
concentric circular boundaries, the effect of using the same control 
functions for the general case will be qualitatively the same, with even 
closer spacing near inner boundary with stronger curvature. Thus the 
control functions derived in the above manner can be expected to produce 
the type of spacing desired in general applications. A version of the 
TOMCAT code incorporating several of these functions has been written 
and has been used to produce coordinate systems for airfoils with the 
spacing at the airfoil set at automatically through (18) . An 

/r 

example is shown In Fig. 4, using the function above (20). Other exam- 
ples are given in Ref. 9. 

D. Revised Generating System 

The form of the control function Q taken in (13) naturally leads 
to the idea of replacing the original elliptic system, (7), with the 
system 

V 2 £ = (C 2 + £ : Z ) p a,n) (35a) 

x y 

v 2 n = (n x 2 + n y 2 ) Q (£,n) (35b) 
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a 


since the terms multiplying P and Q here are, respectively, equal to 
and • With this system the transformed equations are 


ax^ - 23x^ 


yx 


n n 


- (aPx + 


yQ x n ) 


(36a) 


ay^ - 2gy^ + yy^ = - (ctPy^ + yQy^) (36b) 

This form has also been given by Shanks and Thompson, Ref. 10, and by 
Thomas and Middlecoff, Ref. 11. This form has now been adopted in the 
latest version of the TOMCAT code. 

The exponential forms of the functions P and Q, and the discussion 
given therewith above, are still applicable with this system. Appropriate 
values of the attraction amplitudes are several orders of magnitude smaller 
with this new system because of the relatively large values attained by 
the terms multiplying P and Q for small Jacobians. 

Finally, it is useful to solve (36) simultaneously to display P and 
Q explicitly as 


p = “ aj (y n Dx " x n Dy) (37a) 

Q = YJ (y c Dx ~ x , Dy > (37b) 


Dx = ax__ - 

- 2f5x 

+ yx 

(38a) 



nn 

Dy = ay ^ ■ 

■ 26y Cn 

+ yy 

nn 

(38b) 


With (37) the control functions required to produce any specified 
solution x(£,n), y(£ >n) could be calculated. Although such a procedure 
is normally of only academic interest, since the solution x(£,n), y(£,q) 
is yet to be determined, it might be useful in some cases to determine 


P and Q from (37) for some approximate solution generated, say, by 
simple interpolation from the boundaries, and then to use smoothed 
values of these functions as the control functions for the actual 
solution. Although the approximate solution might have lacked continu- 
ity of derivatives, the actual solution determined by solving the 
elliptic system with the smoothed control functions will have continuous 
derivatives, while following generally the form of the approximate 
solut ion. 

E . Control Functions for Near Orthogonality at Boundary 

Another example of the usefulness of (37) is as follows. The 
solution for the concentric circle case can be generalized slightly to 
include variable spacing of points along the boundaries by taking, 
instead of (11) , 


x = r(n) cos 


[ 2 tt 


e (O - 8(1) , 

g(D - gU) 1 


(39a) 


y = r(n) sin 


[2n 


g oh 

g(D - g(D J 


(39b) 


Substitution of these functions in (37) then results in 


P 



(40a) 


Q 



(40b) 


The second of these is the same as (13), using (14) and considering the 
above re-definition of Q, and was used above to generate the control 
function Q. 

With g(n) determined by the boundary point spacing, the control 
function P given here will maintain the £-lines as radial lines, i.e.. 


normal to the circular boundaries. Note that arc length along the 
circular boundary is given by 


8(0 


2tt 


g(C) - b( 1) 
g(I) - g(l) 


r 


( 41 ) 


so that the function g(£) may be related to arc length by 
g(£) = g(l) + ^■ I - 2 ^ r S - ( - 1 -- 8(C) 


(42) 


Thus (40a) can be rewritten in terms of arc length as 



As discussed above for the control function Q, this idea can be 
carried over to the case of general boundaries to produce the same effect 
qualitatively. Thus in the general case, the control function P could 
be determined at each boundary from (43) , and then values of P in the 
field could be taken from linear interpolation between the values at 
corresponding boundary points. 

F. Attraction to Fixed Lines in Physical Space 

As mentioned above, the attraction of coordinate lines to fixed lines 
and/or points in physical space, rather than to floating coordinate lines 
and/or points, requires further consideration. Recall that in the above 
discussion, n-lines are attracted to other p-lines, and £-lines are 
attracted to other £-lines. It is unreasonable, of course, to attempt to 
attract p-lines to £-lines, since that would have the effect of collapsing 
the coordinate system : 



When, however, the attraction is to be to certain fixed lines in 
x-y space, defined by curves y = f (x) , care must be exercised to avoid 
attempting to attract n or £ lines to specified curves that cut the n 
or £ lines at large angles. Thus, in the figure below: 



it is unreasonable to attract £ lines to the curve f (x) , while it is 
natural to attract the n-lines to f (x) . 

However in the general situation, the specified line f (x) will not 
necessarily be aligned with either a £ or n line along its entire length. 
Since it is unreasonable to attract a line parallel to itself, some 
provision is necessary to decrease the attraction to zero as the angle 
between the coordinate line and the given line f (x) goes to zero. This 
can be accomplished by multiplying the attraction function by the cosine 
of the angle between the coordinate line and the line f (x) . It is also 
necessary to change the sign on the attraction function on either side 
of the line f (x) . This can be done by multiplying by the sine of the 
angle between the line f (x) and the vector to the point on coordinate 
line . 

These two purposes can be accomplished as follows. Let a general 
point (x, y) be located by the vector R(x, y) , and let the attraction 
line y = f (x) be specified by the collection of points S(x^, y^) > 
i = 1, 2, — , n. Let the unit tangent to the attraction line be 
t(x i? y^) , and the unit tangent to a £-line be Then the sine 

and cosine of the angle between the £~line and the attraction line may 
be written as 


sine 


[t 


i x (R - S t )] 

" I?'- sj 


k 
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cosine = 


T 

-i 


where k is the unit vector normal to the two dimensional plane. These 
relations are evident from the figure 



The control function P(x, y) may then be logically taken as 


P(x,y) * -T, a (t 


i-1 


in 


T.«>) 

'1 


[ t. x (R - S.)] 

|R - sj 


- X P (~d^ j R - S ± \) 


(44a) 


with the analagous form for Q: 


Q(x,y) = jJ 1 a i (t 1 


. . [t. x (R - S.)] 

(n) -i - ~i 

h } TF^sTi 


-exp (-d^ | R - S.j) 


(44b) 


These functions depend on x and y through both R and or and 

thus must be recalculated at each point as the iterative solution of (36) 
proceeds. This form of coordinate control will therefore be more expres- 
sive than that based on attraction to other coordinate lines. 


There is no real distinction between "line" and M point ,T attraction 
with this type of attraction. ’’Line" attraction here is simply attraction 
to a group of points that form a line f (x) . If line attraction is speci- 
fied, then the tangent to the line f (x) is computed from the adjacent 
points on the line. If point attraction is specified, then the "tangent" 
must be input for each point. 

The tangents to the coordinate lines are computed from 


T^= — (ix + jy ) (45a) 

JZ n ~ n 


X (o) = — (ix_ + jy f ) (45b) 
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G . Point Distribution on Boundary According to Curvature 

One final technique to mention concerns the placement of points 
along the boundary according to the local boundary curvature. Let a 
boundary curve be described by the function y * f (x) . Then if s is arc 
length along the boundary we have 


Now take the rate of change of arc length with the curvilinear 
coordinate, £, along the boundary to be exponentially dependent on the 
local radius of curvature, r, of the boundary. Thus let 


ds 


1 


-br 

- e 


where b is a free parameter. This function causes the arc length to 
change slowly with £ where the curvature is large. 


Then 


dg _ ds ,ds 
dx dx dg 




v- 


assuming x is normalized to vary from 0 to 1 and g to vary from 1 to I. 
The quadrature may be taken numerically if necessary. 

Then for I number of g-points g * 1, 2, — , I on the boundary, the 
corresponding values of x can be determined by inversion of g(x), done 
by interpolation of tabular values if necessary. The arc length between 
each of these points can then be calculated and the value of the free 
parameter b can be adjusted iteratively to produce, say, a specified 
maximum arc spacing along the boundary, or perhaps, to match a specified 
arc spacing at either end or, for that matter, at any given point. In 
application to airfoils, this procedure is applied to the upper surface 
with b chosen to match a specified maximum arc spacing. A separate 
application is then made to the lower surface with b there being chosen 
to match the arc spacing adjacent to the leading edge on the upper 
surface . 

This procedure produces a smooth point distribution on the boundary, 
with points concentrated in regions of large curvature, yet free of the 
rapid spacing changes that lead to coordinate-system-introduced trunca- 
tion error of the type discussed in an earlier section. 
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IV. SOME RECENT APPLICATIONS OF COORDINATE SYSTEMS 


In addition to extensive application to airfoils, as illustrated 
in Fig. 4, in which the transformed plane is an empty rectangle, some 
more general configurations have recently been treated using a trans- 
formed plane that contains rectangular voids as discussed in Ref. 12. 
For example, a coordinate system used in a simulation of a nuclear 
reactor cooling system is shown in Fig. 5, taken from Ref. 13, and 
systems for Charleston harbor (Ref. 14) and a portion of Lake 
Ponchatrain are shown in Figs. 6 and 7. 


V. CONCLUSION 

Control of the spacing of coordinate lines so as to resolve large 
gradients in numerical solution of partial differential equations 
continues to be of paramount importance. Research has provided some 
means of control and of error estimation. The experience gained thus 
far has indicated the versatility of the coordinate systems generated 
from elliptic systems and the possibility of optimization of such 
systems in adaptation to the nature of particular partial differential 
systems and boundary configuration. 
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Figure 4.- Airfoil coordinate system with concentration of lines 
(supplied by Dr. Krishna Devarayalu of the Boeing Company). 
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